#### FIGURE E.1: INDIVIDUAL NEGATIVE CONSEQUENCES
#### Pooled effects for welfare, crime, and employment outcomes

rm(list = ls())
source("./2_code/00_setup.R")

#### STUDY 1 ####

data <- fread(paste0(data_path, "data_study1.csv"), header = TRUE)
data$worked <- as.factor(data$worked)

# Welfare
lin_welfare_economic1 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian1 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic1 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_crime_humanitarian1 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic1 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_employ_humanitarian1 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])


#### STUDY 2 ####

data2 <- fread(paste0(data_path, "data_study2.csv"), header = TRUE)
data2$worked <- as.factor(data2$worked)

# Welfare
lin_welfare_economic2 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian2 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic2 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_crime_humanitarian2 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic2 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_employ_humanitarian2 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])


#### STUDY 3 ####

data3 <- fread(paste0(data_path, "data_study3.csv"), header = TRUE)
data3$employed <- as.factor(data3$employed)

# Welfare
lin_welfare_economic3 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian3 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic3 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_crime_humanitarian3 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic3 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_employ_humanitarian3 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])


#### META ANALYSIS ####

dataset_names <- c("data", "data2", "data3")

### Observation counts for welfare ###

calc_obs_eco_welfare <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_welfare)))))
}

calc_obs_hum_welfare <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_welfare)))))
}

obs_eco_welfare <- unlist(lapply(dataset_names, function(name) calc_obs_eco_welfare(get(name))))
obs_hum_welfare <- unlist(lapply(dataset_names, function(name) calc_obs_hum_welfare(get(name))))


### Observation counts for crime ###

calc_obs_eco_crime <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_crime)))))
}

calc_obs_hum_crime <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_crime)))))
}

obs_eco_crime <- unlist(lapply(dataset_names, function(name) calc_obs_eco_crime(get(name))))
obs_hum_crime <- unlist(lapply(dataset_names, function(name) calc_obs_hum_crime(get(name))))


### Observation counts for employment ###

calc_obs_eco_employ <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_employment)))))
}

calc_obs_hum_employ <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(immigrant_employment)))))
}

obs_eco_employ <- unlist(lapply(dataset_names, function(name) calc_obs_eco_employ(get(name))))
obs_hum_employ <- unlist(lapply(dataset_names, function(name) calc_obs_hum_employ(get(name))))


### Economic treatment, welfare outcome ###

out_economic_welfare <- rbind(
  tidy(lin_welfare_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic_welfare <- cbind(out_economic_welfare, data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco_welfare
))

m.gen_eco_welfare <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_welfare, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_welfare"
)


### Humanitarian treatment, welfare outcome ###

out_humanitarian_welfare <- rbind(
  tidy(lin_welfare_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian_welfare <- cbind(out_humanitarian_welfare, data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum_welfare
))

m.gen_hum_welfare <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_welfare, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_welfare"
)


### Economic treatment, crime outcome ###

out_economic_crime <- rbind(
  tidy(lin_crime_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic_crime <- cbind(out_economic_crime, data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco_crime
))

m.gen_eco_crime <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_crime, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_crime"
)


### Humanitarian treatment, crime outcome ###

out_humanitarian_crime <- rbind(
  tidy(lin_crime_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian_crime <- cbind(out_humanitarian_crime, data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum_crime
))

m.gen_hum_crime <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_crime, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_crime"
)


### Economic treatment, employment outcome ###

out_economic_employ <- rbind(
  tidy(lin_employ_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic_employ <- cbind(out_economic_employ, data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco_employ
))

m.gen_eco_employ <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_employ, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_employment"
)


### Humanitarian treatment, employment outcome ###

out_humanitarian_employ <- rbind(
  tidy(lin_employ_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian_employ <- cbind(out_humanitarian_employ, data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum_employ
))

m.gen_hum_employ <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_employ, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_employment"
)


#### PREPARE PLOT DATA ####

# Helper function to extract forest data from meta object
extract_forest_data <- function(meta_obj, narrative, outcome) {
  total_obs <- sum(meta_obj$n.e)
  data.frame(
    Study = c(meta_obj$studlab, "Pooled Effect"),
    Coefficient = c(meta_obj$TE, meta_obj$TE.random),
    Observations = c(meta_obj$n.e, total_obs),
    SE = c(meta_obj$seTE, meta_obj$seTE.random),
    Lower95 = c(meta_obj$lower, meta_obj$lower.random),
    Upper95 = c(meta_obj$upper, meta_obj$upper.random),
    Lower90 = c(meta_obj$TE - (1.645 * meta_obj$seTE), meta_obj$TE.random - (1.645 * meta_obj$seTE.random)),
    Upper90 = c(meta_obj$TE + (1.645 * meta_obj$seTE), meta_obj$TE.random + (1.645 * meta_obj$seTE.random)),
    Narrative = narrative,
    Outcome = outcome
  )
}

# Extract all forest data
forest_data_eco_welfare <- extract_forest_data(m.gen_eco_welfare, "Economic Narrative", "Welfare Burden")
forest_data_hum_welfare <- extract_forest_data(m.gen_hum_welfare, "Humanitarian Narrative", "Welfare Burden")
forest_data_eco_crime <- extract_forest_data(m.gen_eco_crime, "Economic Narrative", "Crime")
forest_data_hum_crime <- extract_forest_data(m.gen_hum_crime, "Humanitarian Narrative", "Crime")
forest_data_eco_employ <- extract_forest_data(m.gen_eco_employ, "Economic Narrative", "Employment")
forest_data_hum_employ <- extract_forest_data(m.gen_hum_employ, "Humanitarian Narrative", "Employment")


#### PLOTTING ####

desired_order <- c("Pooled Effect", "Study 3", "Study 2", "Study 1")

# Helper function to create forest plot
create_forest_plot <- function(forest_data, caption_text) {
  forest_data <- forest_data %>%
    mutate(Study = factor(Study, levels = desired_order))
  
  ggplot(forest_data, aes(x = Coefficient, y = Study)) +
    geom_point(aes(color = Study), size = 3) +
    geom_errorbarh(aes(xmin = Lower95, xmax = Upper95, color = Study), height = 0, linewidth = 0.9) +
    geom_errorbarh(aes(xmin = Lower90, xmax = Upper90, color = Study), height = 0, linewidth = 1.9) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    labs(title = "", x = "", y = "", caption = caption_text) +
    theme_minimal() +
    theme(
      axis.title.x = element_blank(),
      plot.caption = element_text(size = 11, hjust = 0.5),
      strip.text = element_text(size = 12),
      legend.position = "none"
    ) +
    scale_x_continuous(limits = c(-0.3, 0.12)) +
    scale_color_manual(values = c("Study 1" = "grey50", "Study 2" = "grey50", "Study 3" = "grey50", "Pooled Effect" = "red")) +
    facet_grid(~Narrative)
}

# Create individual plots
forest_welfare <- rbind(forest_data_eco_welfare, forest_data_hum_welfare)
forest_crime <- rbind(forest_data_eco_crime, forest_data_hum_crime)
forest_employ <- rbind(forest_data_eco_employ, forest_data_hum_employ)

coefficient_plots_welfare <- create_forest_plot(forest_welfare, "ATE Welfare Burden")
coefficient_plots_crime <- create_forest_plot(forest_crime, "ATE Crime")
coefficient_plots_employ <- create_forest_plot(forest_employ, "ATE Employment")

# Combine all plots
combined_plot <- coefficient_plots_welfare / coefficient_plots_crime / coefficient_plots_employ

combined_plot

ggsave(filename = paste0(plot_path, "figure_E1.png"), plot = combined_plot, width = 7, height = 9, dpi = 600)

